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We propose and study a simplified model for the surface and bulk structures of crosslinked polymer 
gels, into which voids are introduced through templating by surfactant micelles . Such systems wer e 



recently studied by Atomic Force Microscopy [M. Chakrapani et al, e-print cond- mat/0112255 
The gel is represented by a frustrated, triangular network of nodes connected by springs of ran- 
dom equilibrium lengths. The nodes represent crosslinkers, and the springs correspond to polymer 
chains. The boundaries are fixed at the bottom, free at the top, and periodic in the lateral direction. 
Voids are introduced by deleting a proportion of the nodes and their associated springs. The model 
is numerically relaxed to a representative local energy minimum, resulting in an inhomogeneous, 
"clumpy" bulk structure. The free top surface is defined at evenly spaced points in the lateral (x) 
direction by the height of the topmost spring, measured from the bottom layer, h{x). Its scal- 
ing properties are studied by calculating the root-mean-square surface width and the generalized 
increment correlation functions Cq{x) = {\h{xo -I- a;) — h{xo)\'')- The surface is found to have a 
nontrivial scaling behavior on small length scales, with a crossover to scale-independent behavior on 
large scales. As the vacancy concentration approaches the site-percolation limit, both the crossover 
length and the saturation value of the surface width diverge in a manner that appears to be propor- 
tional to the bulk connectivity length. This suggests that a percolation transition in the bulk also 
drives a similar divergence observed in surfactant templated polyacrylamide gels at high surfactant 
concentrations. 

PACS numbers: 82.70.Gg, 61.43.Hv, 68.37.Ps, 89. 75. Da 



I. INTRODUCTION 



CrossUnked polymer hydrogels are completely interconnected polymer networks that combine high water content 
with high porosity, forming macroscopic molecules that have applications in many fields. Their wide range of pore 
sizes makes such gels ideal for separation of biological macromolecules by electrophoresis or chromatography . The 
crosslinkingprocess induces a reorganization of the polymer structure, resulting in inhomogeneities in the spatial 
density l^j, ^ |[ ||, |^. These inhomogeneities affect the surface configurations of the gels, but only in recent years 
have advances in atomic force microscopy (AFM) made possible the imaging of soft material surfaces in an aqueous 
environment 0|. The effects of the crosslinking density, temperature, pressure, and sample thickness on the surface 
topography have been studied by Suzuki et al. |^ • Their results indicate that the structural features of the surface 
on both the micrometer and nanometer scales depend on these factors. It has been suggested Q that control of the 
characteristic length scale of the gel surface using external stimuli may have applications in a variety of fields, such 
as regulation of adsorption and release of specific molecules by the intermolecular forces between the surface and 
the molecule. Thus, gel surfaces provide a unique opportunity to explore the interplay between phenomena on the 
macroscopic and nanoscopic scales. 

Templated polyacrylamide gels are formed by polymerizing acrylamide with a crosslinker in the presence of a 
surfactant. The surfactant molecules form monodisperse micelles of a size roughly comparable with the crosslinker 
separation, and the presence of these micelles alters the gel pore structure, enhancing the gel's separation properties 
. In a recent work the surface morph olog y of templated polyacrylamide gels was extensively studied by AFM and 
scaling analysis of the resulting images . This study indicates that the gel surfaces are self-affine on short length 
scales, with roughness (Hurst) exponents on the order of 0.8-1. In the absence of surfactant a cross-over length, above 
which the surface is no longer self-affine, was estimated to be on the order of 300 nm, and the saturation value of the 
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interface width was on the order of 1 nm. Both values increased dramatically with the introduction of surfactant. 

Detailed kinetic lattice models of the polymerization of crosslinked polymer gels have previously been constructed 
pTf . However, they emphasize the kinetics of polymerization and are extremely computationally intensive. Inspired 
by the experiments reported in Ref. |l0[] , the aim of the present work is rather to construct and study a simple 
continuum model that can reproduce some of the observed scaling characteristics of the templated gel surfaces at a 
more modest (but still substantial) computational cost by concentrating on the elastic structure of the gel. Some 
preliminary results were presented in Ref. 

The remainder of this paper is organized as follows. In Sec. II we introduce the model and detail some aspects of 
the numerical calculations. In Sec. Ill we recall some scaling concepts associated with the analysis of surfaces. In 
Sec. IV we show the results of our calculations, in particular the scaling properties of the surface width and increment 
correlation function along with their dependence on the size of the sample. We also comment on the relevance of our 
results for the interpretation of recent AFM experiments. Finally, in Sec. ^ we present our conclusions. 



II. MODEL 



The model consists of a two-dimensional network of nodes interconnected by massless springs. The nodes represent 
crosslinker molecules, and the connecting springs represent polymer chains. The network topology consists of a 
triangular lattice of nodes, each of which is connected by harmonic springs to its six nearest neighbors (except at the 
top and bottom surfaces, where each node has only four connections). The model has no excluded- volume interactions. 
A triangular lattice was chosen to ensure geometrical stability in two dimensions without introducing bond-angle 
constraints. The corresponding unrealistically high crosslinker functionality should not significantly influence our 
results. The network has periodic boundary conditions in the horizontal (x) direction, free boundary conditions at 
the top {y > 0) layer, and the nodes in the bottom layer are fixed at y = 0, corresponding to bonding to a rigid 
substrate. There are Lx nodes in the horizontal direction and Ly nodes in the vertical direction. 

The total energy of the network is E = X]i(l/2)^i(^i ~ 'oi)^) where ki, li, and loi, are the spring constant, the actual 
length, and the equilibrium length under zero external force of the ith spring, respectively. (All quantities in this 
paper are given in dimensionless units.) The equilibrium length of each spring is independent of the other springs 
and is randomly chosen with probability density function (pdf) 

P{lM)^2-flo^exp{~-fll^) , (1) 

where 7 is proportional to the inverse of the average number of monomers between crosslinkers. This pdf corresponds 
to the case that the equilibrium distance between crosslinkers is proportional to the square root of an exponentially 
distributed number of monomers. It is consistent with the picture that crosslinkers are distributed randomly along the 
polymer chains, and that a spring of equilibrium length Igi corresponds to a polymer of the same average end-to-end 
distance in the random-coil collapsed phase p3| , p^ . The average equilibrium length of a spring is (/oi) = \/7iy7/2. In 
agreement with a mean-field spin-chain approximation for the elastic properties of a collapsed polymer chain p4| , we 
require that the elastic constant of the ith spring should be inversely proportional to its equilibrium length, ki — l^^ 
in our dimensionless units. 

With the above distribution of equilibrium lengths, the pdf for the force F exerted by a spring of length I is 
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(2) 



1 -F 

for < 1, while it vanishes for F > 1. The average force exerted by a spring of length / is then 

{F\l) (3) 

Initially, springs with equilibrium lengths drawn independently according to P{lo) are placed on the bonds of a 
regular triangular lattice of unit lattice constant, such that all of them are stretched or compressed to a length of 
one. In order that this initial configuration should not be globally stressed, we impose the condition that the average 
force exerted by a spring of unit length should be zero. By Eq. (^), this is satisfied for 7 = I/tt. However, the initial 
configuration is locally stressed, since h — 1 for all i while (/oi) = 7r/2 « 1.57. The system is frustrated in the sense 
that there is no configuration that can simultaneously minimize the energies of all the springs. As a result there is a 
large number of configurations corresponding to local energy minima with similar energies. 

The network was relaxed using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) quasi-Newton 



minimization algorithm |15, Q6| until a locally stable configuration was reached, corresponding to a local minimum 



in the global energy landscape. By minimizing several statistically equivalent realizations of the random spring 
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configuration, we checked that the local minima reached by this procedure are representative of typical relaxed 
configurations. Although the total energies of different relaxed configurations were not completely identical, all the 
structural results were the same to within our numerical accuracy. 

The L-BFGS subroutine requires the calculation of the total energy and its gradients, i.e., the local forces in the 
vertical and horizontal directions for each node. By monitoring the total energy and the absolute value of the force we 
found that the subroutine must be called about 3 x 10^ times to ensure that the system reaches a local minimum. Since 
each relaxation step requires a new evaluation of the total energy and the local forces, the L-BFGS algorithm is quite 
computationally intensive, taking about 1-2 days on a current 600 MHz dual-processor Pentium III PC with 512 MB 
of shared memory to relax a typical system with about 10^ degrees of freedom. However, it was far more efficient for 
our application than a naive steepest-descent method, requiring about an order of magnitude less computer time while 
leading to essentially equivalent final energies and structures. After relaxation, we examined the bulk and surface 
properties of the network. Images of a typical relaxed network are shown in Figs. |l| and |[ Inhomogeneities in the 
spatial crosslinker density are evident in these figures. 



III. SCALING PROPERTIES OF SURFACES 



The scaling properties associated with the height fiuctuations of rough surfaces have been the object of a great deal 
of interest. Vapor deposition, epitaxial growth, polymer growth, and many other phenomena exhibit scale invariance 
[T^ , [l9| , po[ . In a statistical sense, the height of a single- valued self-affine random surface, h{x), satisfies the scaling 
relation 

h{x) = X-"h{Xx) , (4) 

where A is a dimensionless scaling parameter and H is the so-called Hurst exponent, which in this case is identical 
to the roughness exponent pof . Consequently, the increment correlation function C{x) of a translationally invariant, 
self-affine surface scales as 

C{x)^([h{xo + x)-h{xo)]^)^x^" , (5) 
where the average is taken over Xq. Many surfaces present a mixed behavior characterized by 

C(x) ~ I ^'^^ ^ (6) 

^ ' 1 constant for a; >> ' 

where ly, is a cross-over length scale separating the two behaviors. However, there exist surfaces for which the 
scaling cannot be described by a single exponent; instead an infinite hierarchy of characteristic exponents has to 
be introduced. These surfaces are called multi-affine, and they are occasionally observed in nature ||2l|, p2| , p3| . 
Molecular-beam epitaxy (MBE) models ||2^, polished metal surfaces geomorphological features |26|, and time 



series of financial data |27, Eq| provide some examples of multi-affine scaling. The characteristic exponents of a 



multi-affine function can be calculated from the gth-order increment correlation function, defined as 

Cq{x) = {\h{x^ ^x)- /i(xo)r) ^ (7) 

with the generalized Hurst exponent Hq depending continuously on g, at least in some region of q values. Thus C2{x) 
is identical to the usual increment correlation function C{x) defined in Eq. (|^). Like self-affine surfaces, multi-affine 
surfaces can show a cross-over to scale-independent behavior, analogous to Eq. (||). 

In the next section we present our results on the structure of the bulk and the scaling behavior of the surface, 
calculated after the spring network has been relaxed to a representative local energy minimum. 



IV. NUMERICAL RESULTS 



We analyzed spring networks of sizes Lx x -^y In order to simulate experimental measurements of the scaling 
behavior of surfaces of gels in which pores have been introduced by polymerizing acrylamide plus a crosslinker 
templated with surfactants that are later removed we created voids in the networks. The voids were introduced 

by selecting a proportion of the nodes at random and removing these nodes and the springs associated with them. 



Next, the Hoshen-Kopelman algorithm |29, 30 was applied to identify and discard clusters that were not connected 



to the fixed substrate. This last step is particularly important near the site-percolation limit, which for the triangular 
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lattice is at a vacancy concentration of 50% Q. Note that by "vacancy concentration" we always mean the 
initial concentration of vacant nodes before removal of the disconnected clusters. The total volume fraction of the 
void structure resulting from the cluster removal is of course larger. 

We define the onc-dimcnsional surface as the set of surface heights, h{xj), at equally spaced discrete points, Xj G 
[0,Lx), considering the polymer chains as straight line segments between nearest-neighbor connected nodes. To 
calculate the surface height at Xj, we first identify the set of all springs which intersect a vertical line at Xj. From 
this set of springs, we define 

h{xj) =ma-x[yi{xj)], (8) 

where yi{xj) is the vertical location of the intersection of the ith spring at xj. We use the sign convention that y = 
for the fixed bottom nodes and the free surface is located at positive values of In Fig. |l| we show details of the bulk 
and free surface of systems with vacancy concentrations of 0%, 40%, and 49.5%, while in Fig. || we show a magnified 
view of a section of the bulk of a vacancy-free system. 

To obtain an estimate of the finite-size effects and select an appropriate system size for the production runs, we 
first considered vacancy-free spring networks of sizes x Ly, with = 512, 1024, and 2048, and Ly = 256, 384, 
512, 640, and 748. We calculated the root-mean-square width of the surface (also known as the rms roughness) on a 
length scale L, wl, by recursively subdividing the surface into n nonoverlapping segments of length L = Lx/n, with 
n an integer between 2 and Lx/2. The mean-square width of the fcth segment is 

wlk - (y'>Lfc - (y)lk , (9) 

where {■)Lk denotes averaging within the segment. The rms width is defined as the square root of the average of w^j, 
over the n segments, vol = Z]fc=i "^Lk- ^^'^ statistical errors in wl were estimated in the usual way from the 

standard deviation of w\^.. 

At this point we emphasize that the order in which the average and the square root are taken in the definition 
oi Wl is important. For a simple self-affine surface it is irrelevant whether the square root is taken before or after 
the average; in either case one obtains the same exponent. However, for a multi-affine surface the results are quite 
different. In this case the exponent obtained from wl is equal to that calculated from the increment correlation 
function C(V) only when the average is calculated before taking the square root, as it is done here. The results for 
Wl, averaged over four independent configurations for each system size, are shown in Fig. |^. 

The data shown in Fig. ^ indicate that, for the values of and Ly considered, wl is independent of the system 
size. However, when comparing the average height for the jth layer in the relaxed system, {yj), with its original 
value in the regular triangular lattice, we found a slight contraction limited to approximately the top 250 layers. To 
minimize finite-size effects due to interference of the finite thickness Ly with this skin effect, while maintaining L^ as 
large as possible so as to reach large horizontal length scales and reduce statistical uncertainties in the estimates of 
wl, we chose a system size of L^ x Ly = 1024 x 768 for the production runs. Both the bulk and surface structures 
were analyzed as discussed in the following two subsections. 



A. Bulk structure 



In this subsection, we discuss the bulk properties of the spring network far from both the fixed and free surfaces. 
We only present results for a gel without vacancies. For gels with high vacancy concentrations, it is difficult to define 
a bulk region, since the width of the surface region increases with increasing vacancy concentration. To study the 
bulk structure in this case, simulations with periodic boundary conditions would be preferable. 

The presence of inhomogeneities in the spatial crosslinker distribution, where the crosslinkers and the polymer 
chains (shown as line segments) have greater densities in some regions (clumps), can be seen in Fig. g and Fig. |. 
The effect is more pronounced for higher vacancy concentrations. This is a desirable feature of the model since static 
inhomogeneities are characteristic of crosslinked polyacrylamide gels and have been related to a freezing-in of the 
topological structure during polymerization |3[ ||, ||] . For gels without vacancies. Fig. |l](a), the clumping is not as 
obvious to the eye. Thus, we show in Fig. |^ a magnified image of a small region of the bulk of a vacancy-free system, 
in which the clumping can be clearly seen. To quantify the clumpirig, we characterize the bulk scaling properties of 
the gel by a method analogous to the width-scaling method of Eq. (y) . 

We define p{x, y) — 5(x — Xi,y — yi) to be the density of nodes at a particular point {x, y) of the bulk. Here 
6{-) is the two-dimensional Dirac delta function, and the sum runs over all nodes {xi,yi) in the system. We calculate 
the mean and the variance of the number of nodes inside a box of side L, {N)l and {N'^)l — respectively, 
by averaging p{x, y) over non-overlapping boxes. When there are no correlations within the gel, corresponding to a 
purely random distribution of node positions, the number of nodes in a box obeys a Poisson distribution. This yields 
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{N)l/{{N'^)l — = 1. However, when the nodes have a more regular distribution, the denominator is reduced 

and the ratio increases. In Fig. || {N)l/{{N'^)l — (^)i) — 1, averaged over eight independent reaUzations of the 
system, is shown vs L from L = 0.25 to 10. For small L, the function tends to zero, indicating that nodes which are 
very close together are essentially randomly placed. The function goes through a maximum near the average spring 
length (L « 1.57), corresponding to a more regular distribution of the nodes on this scale. 



B. Surface scaling behavior 

We calculated the qih root of the gth-order increment correlation function, Cy'^(L), for q = 0.5, 1, 1.5, 2, and 2.5. 
The results for systems with 0%, 40%, and 49.5% vacancies are shown in Fig. |^(a), (b), and (c), respectively. For each 
of the two lower vacancy concentrations the data were averaged over ten independent realizations of the system, while 
thirteen realizations were used for the highest concentration. For all three vacancy concentrations, the surface is seen 
to be multi-affine with q-dependent slopes for L below a crossover length, ly . [For the highest vacancy concentration 

(Fig. |^(c)), two scaling regimes with different, nonzero Hq are seen for q < 1.] For L > Ix, Cg^' reaches a saturation 

value Cl^it ■ Both of these length scales, /x and C^l^ , depend on the vacancy concentration and q in ways that are 
discussed below. 

By numerically calculating the logarithmic derivatives of the correlation functions with respect to L we determine 
the scaling exponents, 

q a log i 

Effective values for Hq for different L and q were estimated as two-point derivatives over intervals corresponding to 
a doubling of L. The results are shown in Fig. ||(a), (b), and (c) for systems with 0%, 40%, and 49.5% vacancies, 
respectively. The figures clearly indicate that on length scales below Ix , the surfaces are multi-affine with g-dependent 
scaling exponents Hq. In this scaling regime the effective exponents are approximately independent of L. Scaling 
exponents Hq for different values of the vacancy concentration are presented in Table |. The exponents shown in 
this table were obtained from Eq. (^0|) at fixed L = 0.088, in the small-L scaling regime. In Fig. ^ we show how 
these values of Hq change with q and the vacancy concentration. For g > 1, Hq decreases with increasing vacancy 
concentration, while it increases for g < 1. For q = I, Hq is slightly below unity and independent of the vacancy 
concentration to within our numerical accuracy. 

Multi-affinity is relatively rare in natural surfaces, and, in fact, it is not observed in the experimental system that 
inspired the present study ||l0| . It is therefore reasonable to inquire about the source of the multi-affinity observed in 
our model system. We find that it is caused by the finite density of discontinuities, which result from overhangs in 
the simulated surfaces (see Fig. 0) . This effect is illustrated in two ways in Fig. ^. In Fig. ||(a) we show on a log- log 
scale Cg(L)^/* for different values of q for surfaces generated by fitting linear line segments of unit length to self-affine 
surfaces with a Hurst exponent of 0.75, which were generated by the method of successive additions |32j. These 
surfaces are seen to be self-affine with a Hurst exponent which crosses over from near unity to approximately 0.75 at 
about L = 1. In the same figure we also show (by thin curves) the corresponding generalized increment correlation 
functions for surfaces generated from the above ones by the insertion of a finite density of vertical discontinuities. 
In contrast to the result for the original surfaces, the slopes of the latter curves depend on q, indicating that the 
discontinuous surfaces are multi-affine. Conversely, in Fig. |g(b) we illustrate the effect of removing the discontinuities 
from a surface generated by our model with 0% vacancies. While the original surface is multi-affine (signified by the 
g-dependent slopes), the surface from which the discontinuities are removed is self-affine with a Hurst exponent of 
approximately 0.98, corresponding to the straight line segments used to represent the polymer chains. These numerical 
results clearly show that the multi-affinity in our model surfaces is due to the vertical discontinuities caused by the 
overhangs. Possible methods to improve the agreement of the scaling behavior of the model surfaces with that of the 
experimental system by removing the discontinuities are discussed in Sec. ^ 

The large-L saturation values, , are summarized in Fig. ^. For a fixed concentration of vacancies, the satu- 
ration value increases with increasing q. The saturation value also increases with the vacancy concentration, quite 
dramatically so when the concentration approaches the site-percolation limit of 50%. Figure ^(a) shows C^l^ vs the 
vacancy concentration on a linear-log scale, while Fig. ^(b) shows the same quantity vs the relative distance from 
the percolation limit, e = 1 — (%vacancies)/50%, on a log-log scale. The straight line in Fig. ^(b) is proportional to 

e~'^ with V = A/3 being the exact value for the connectivity- length exponent |3l|. The figure indicates that 
appears to be coupled to the connectivity length for vacancy percolation in the bulk, at least for q> 1. 

For a given vacancy concentration, all the curves in Fig. ^ coincide for L ^ l^, indicating that the surface structure 
is not scale-dependent on these large length scales. Due to the difficulty in obtaining a precise estimate for Zx, we 
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chose to define it as the value of L at which Hq = 0.3. This value is shown vs the vacancy concentration in Fig. |lO| 
which, like Fig. ^, is divided into two parts. In Fig. p^a) /x is shown vs the vacancy concentration on a linear- log 
scale, while in Fig. |l^(b) it is shown vs e on a log-log scale. The increase in the vicinity of the percolation threshold 

is similar to that of Cj/t' in Fig. H, except that Ix has a significant background value which dominates at and below 
a vacancy concentration of 40% (e > 0.2). Thus our data indicate that both these characteristic surface lengths 
are coupled to the bulk connectivity length as percolation is approached. Similar behavior to that observed for the 
saturation value of C2{Ly^^ is also seen for the saturation value of wl- 

In the templated polyacrylamide gels studied in Ref. a strong increase in the crossover length and surface 
width at surfactant concentrations above 20% by weight coincides with a change from an optically clear material to a 
white, opaque one. Such a change in the light scattering intensity is consistent with an increase in the bulk correlation 
length from much smaller than, to on the order of the wavelength of visible light. Considering the numerical and 
experimental evidence together, we conclude that the dramatic increases in the surface-related length scales that occur 
with increasing volume fraction of voids are related to a percolation transition in the system of vacancies in the bulk, 
both in our model and in the experimental system by which it was inspired. 



V. SUMMARY AND CONCLUSIONS 



In this paper we have proposed and studied a simple frustrated spring-network model which, with a relatively modest 
computational effort, reproduces aspects of the scaling behavior of the surfaces of surfactant-templated polyacrylamide 
gels, recently observed in AFM experiments These results suggest that simple models that incorporate elastic 
properties of the gel can provide a useful framework for understanding the structural characteristics of these materials, 
which are otherwise extremely difficult to simulate due to their complexity. 

The main similarities and differences between the surface structures of the model and the experimental system are 
as follows. The main similarity is a nontrivial scaling behavior corresponding to a power-law form of the increment 
correlation function at small and intermediate length scales. This power-law behavior terminates at a crossover 
length scale, above which the rms surface width and the increment correlation functions both reach saturation values 
that appear to be independent of the system size. The origin of this crossover behavior is explained in terms of 
the structural inhomogeneities of the gel. On large length scales the scaling behavior reflects the average network 
structure, while on small length scales the microscopic density fluctuations of the gel network play the most important 
role. The quantitative dependence of the crossover on the bulk vacancy concentration is discussed further below. The 
main difference between the model and the experimental surfaces lies in the details of the scaling behavior at length 
scales below the crossover length: the model surfaces are multi-affine, while the experimental surfaces are self-affine. 
The multi-affinity in the model is caused by the finite density of vertical discontinuities in the simulated surfaces 
(see Fig. |l|), which are not seen in the experimental system (see Fig. 7 of Ref. The discontinuities in the 

simulated surfaces as they are defined here result from overhangs in the spring configuration at the surface, and 
modification of the model to improve the detailed agreement with the scaling behavior of the experimental system 
should aim at eliminating these overhangs. (It is also possible that overhangs do in fact exist in the undisturbed 
experimental surfaces, but are so mechanically weak that they are smoothed out by the force of the AFM tip.) In 
the model, the number and size of overhangs could be reduced or eliminated in several ways. One possibility is to 
simulate the smoothing action of the AFM tip, e.g., by a running-average method. Another possibility is to include 
a weak gravitational field (the specific gravity of the polymers is greater than that of water), which would bend 
overhanging sections down into the main part of the interface. It should also be noted that the larger coordination 
number in three dimensions most likely would reduce the importance of overhangs in a three-dimensional version 
of the model, even without further modifications. (Extension to three dimensions would also reduce the number of 
disconnected pieces, which are partly responsible for the large overhangs seen in the case of high vacancy densities. 
However, the computational cost would be large.) A more exhaustive numerical and analytical study of the infiuence 
of discontinuities on surface scaling behavior is currently underway |Q . 

The experimental surfactant templating was simulated by introducing a nonzero volume fraction of voids in the 
network through deleting lattice nodes and their associated springs. Both the horizontal crossover length and the 
vertical surface thickness (as measured by the saturation values of both the rms surface width and the generalized 
increment correlation functions) increase dramatically as the vacancy concentration approaches the site-percolation 
limit of 50%. An analogous increase of the characteristic horizontal and vertical length scales with surfactant concen- 
tration is observed in the experimental system, where it coincides with a dramatic increase in the optical opacity. The 
numerical results for our model thus support the hypothesis that this increase in the experimental system is caused 
by a percolation transition of the voids created by removal of the surfactants. 

In conclusion, despite the differences in the details of the scaling behavior at small and intermediate length scales, 
our simplified random-spring model successfully mimics several of the surface characteristics of an important new 
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class of porous media, including the dramatic increase in surface roughness that occurs with increasing surfactant 
concentration. Efforts to further improve the agreement between the small-scale scaling behavior of our theoretical 
model and that of the experimental system by which it was inspired are left for future study. These efforts should 
in particular involve investigations of the dependence of the scaling properties of the model on external forces due 
to gravity or the AFM tip, as well as on the statistics of the local spring distribution and the dimensionality of the 
model. 
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TABLE I: Scaling exponent Hq : dependence on q and vacancy concentration. 



Vacancies 




til 


^1.5 




rr 


n 0/ 
U 70 


1.0008 


0.9582 


0.8683 


0.7030 


0.5428 


10 70 


1.0113 


0.9649 


U.OZOZ 


0.d2d8 


0.4688 


20 % 


1.0286 


0.9543 


0.7657 


0.5627 


0.42577 


Qn 0/ 
oU /o 


1 nc; /1 1 


O.yo4i 


U. f ziz 


U.OZOD 


0.4UOO 


40 % 


1.0910 


0.9607 


0.7045 


0.5168 


0.4059 


45 % 


1.1087 


0.9588 


0.6940 


0.5105 


0.4019 


47% 


1.158 


0.9713 


0.6939 


0.5105 


0.4019 


49 % 


1.140 


0.9662 


0.6820 


0.5052 


0.4021 


49.5 % 


1.1552 


0.9711 


0.6772 


0.5025 


0.4008 



8 




FIG. 1: Details of the surface and bulk structures of relaxed 1024 x 768 systems with vacancy concentrations of (a) 0%, (b) 
40%, and (c) 49.5%. The thin line segments represent polymer chains. The bold lines (vertically displaced for easier viewing) 
represent the free surface. 
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FIG. 2: Magnified detail of the relaxed bulk structure for a 1024 x 768 system without vacancies. 
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FIG. 3: The rms surface width wl for vacancy-free spring networks, shown on a log-log scale as a function of the length scale 
L for a number of system sizes. The two scaling regimes defined in Eq. are clearly seen. Averaged over four independent 
configurations for each system size. 
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FIG. 4: The bulk scale parameter vs the length scale L. The data are for a system of 1024 x 768 nodes with 0% vacancies, 
averaged over eight independent configurations. The function values near zero indicate randomly placed nodes at small length 
scales, and the maximum at L r; 1.57, which is the average equilibrium spring length, indicates that nodes are most regularly 
placed on this scale. 
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FIG. 5: The gth root of the generalized increment correlation function Cq(L) for systems with vacancy concentrations (a) 0% 
vacancies, averaged over ten independent realizations, (b) 40% vacancies, averaged over ten independent realizations, (c) 49.5% 
vacancies, averaged over thirteen independent realizations. 
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FIG. 6: The effective generaiized Hurst exponent Hq for systems witli different vacancy concentrations, obtained as iogaritlimic 
two-point derivatives of C\^''{L) with respect to L. The data for C\^''{L) were averaged as in Fig. |5| For small L and g, the 
derivative approaches unity, indicating Hq ~ 1.0. For large L, Hq approaches zero, signaling the scale-independent saturation 
regime, (a) 0% vacancies, (b) 40% vacancies, (c) 49.5% vacancies. 
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FIG. 7: The generalized Hurst exponent Hq for systems witii different vacancy concentrations, from Table ^ (a) Shown vs q 
for different vacancy concentrations, (b) Shown vs vacancy concentration for different values of q. In both parts the lines are 
merely guides to the eye. 
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FIG. 8: Effects of discontinuities on the generalized correlation functions Cq{LY^'^ and their logarithmic slopes, the generalized 
Hurst exponents Hq, for self-affine and multi-afhne surfaces. In all sets of curves shown, q goes from 0.5 to 4.0 in increments 
of 0.5, from below to above in the figures. For clarity of view, the data corresponding to continuous surfaces (the thick curves) 
have been divided by ten. (a) The thick curves correspond to continuous, approximately self-affine surfaces, as indicated 
by the q- independent Hq. The thin curves represent discontinuous surfaces obtained by inserting a finite density of vertical 
discontinuities in the original surfaces. For the latter surfaces, Hq depends on q, indicating multi-affinity on length scales less 
than ten. Each curve is an average over data for 1000 surfaces, (b) The thin curves correspond to a multi-affine simulated 
gel surface with 0% vacancies, corresponding to those for which Cq{L)^^'' are shown in Fig. ^a). The thick curves correspond 
to a continuous surface generated from the simulated gel surface by removing all vertical discontinuities. In this case, Hq is 
independent of q, indicating simple self-affinity. 
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FIG. 9: The large-L saturation value of Cq^''(L), . (a) Shown vs vacancy concentration on a hnear-log scale, (b) Shown 
vs e = 1 — (%vacancies)/50% on a log-log scale. (Note that the vacancy concentration here increases toward the left.) The 
straight line is proportional to e"^''^, corresponding to the connectivity length for bulk site percolation In both parts 

the error bars are smaller than the symbol size. 
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FIG. 10: The crossover length Ix between the multi-afhne regime and the saturation regime, estimated as the minimum value 
of L for which Hq{L) = 0.3. (a) Shown vs vacancy concentration on a linear-log scale, (b) Shown vs e = 1 — (%vacancies)/50% 
on a log-log scale. (Note that the vacancy concentration here increases toward the left.) The straight line is proportional to 
e"*'''', corresponding to the connectivity length for bulk site percolation [ ^o| , |3l| ]. In both parts the error bars are smaller than 
the symbol size. 
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